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Abstract 

We analyse the dynamics of expectation values of quantum observables 
for the time-dependent semiclassical Schrodinger equation. To benefit 
from the positivity of Husimi functions, we switch between observables 
obtained from Weyl and Anti-Wick quantization. We develop and prove 
a second order Egorov type propagation theorem with Husimi functions 
by establishing transition and commutator rules for Weyl and Anti-Wick 
operators. We provide a discretized version of our theorem and present 
numerical experiments for Schrodinger equations in dimensions two and 
six that validate our results. 

1 Introduction 

Our investigation is devoted to the time-dependent Schrodinger equation 

iedtiP = -'^^il^ + Vi>, V(0)=^o (1) 

for general square integrable initial data iI^q : ^ C with HV'oIIl^ = 1- We 
assume that the potential 1/ : R'' — >■ R is smooth and satisfies suitable growth 
conditions at infinity, such that the semi-classical Schrodinger operator 

is essentially selfadjoint. Such potentials are met in molecular quantum dy- 
namics for the effective description of nuclear motion after performing the time- 
dependent Born-Oppenheimer approximation |ST01[ Theorem 1]. In this appli- 
cation, the parameter e is the square root of the ratio of electronic versus average 
nuclear mass. Hence, e > is small and typically ranges between i/iooo and i/io. 
This implies that the solution of ([ij is highly oscillatory in space and time with 
frequencies of the order i/e. However, the time evolution of expectation values 

t ^ {lpt,A'ipt)L'^ 

is less oscillatory, and the Egorov theorem |BR021 Theorem 1.2] suggests a 
semiclassical approximation using the Hamiltonian flow $* : M^'' — > M^'* of the 
classical equations of motion q = p, p = — Vy(g). 



For this approximation we adopt the phase space point of view and consider 
operators A = op^®(a), which are obtained from functions a : M^'^ — ;> K by 
the Weyl quantization. Expectation values of Weyl quantized operators can be 
written as phase space integrals 

(V't,opW°(a)^4)i2 = / a{z)W'{iJt){z)dz, 

where W^('!/'t) : K^'^ — )■ M is the Wigner function of the square integrable func- 
tion ■i/'t : R'' — >■ C. The Egorov theorem can be expressed in terms of Wigner 
functions as 

(V't,opW-(a)Vt)L2 - / {ao^'){z)W%^Po){z)dz + 0{e^). (2) 

In general, Wigner functions attain negative values. However, a proper 
smoothing by a Gaussian function with mean zero and covariance | Id 

results in a nonnegative function 

= W^(^t)*G,/2, 

the so-called Husimi function (GMMP97[ Remark 1.4]. The gain in positivity 
is accompanied by a loss in accuracy with respect to time propagation, since 
replacing the Wigner function >V^(?/'o) by the Husimi function ■H'^(V'o) in the 
Egorov theorem (|2| deteriorates the approximation error from second to first 
order in e. This motivates the question, whether an e-corrected flow can com- 
pensate this accuracy loss when propagating Husimi functions. 

Our main result answers the question positively. Its basic constituents are 
the e-dependent modification 

= a — |Aa 

of a phase space function a : M^'* — > M by subtracting its Laplacian, the flow 
: R'^'^ — > M^"* associated with the Hamilton function 

together with the solutions A* e M^''^^'^ and F* e M^'' of two first order ordinary 
differential equations, whose right hand side is built from evaluations of D^h 
and D^h along $*, respectively. We obtain: 

Theorem. There exists a constant C = C{a,V,t) > such that for all square 
integrable initial data ipo : M'^ — > C with HV'oIIl^ = 1 and all e > the solution 
of the Schrodinger equation satisfies 

{tPt,op'^%a)A)L2- f Fl{a,){z)W{i^o){z)dz 

with 

Flia,) = a, o - I (F* • (Va o $*) + tr (A* {D^a o $*))) . 



2 



The method of proof rehes on a systematic transition between Weyl and 
Anti-Wick quantization combined with an e-expansion of commutators in Weyl 
quantization. Therefore, the generalization of the approximation to errors of 
the order e'^, fc > 2, is possible. 

The theorem allows an algorithmic realization in the spirit of the particle 
method developed in |LR10j , which is suitable for the fast computation of ex- 
pectation values in high dimensions: sample the initial Husimi function 'H^ijjjo), 
discretize the ordinary differential equations for $*, and Af, propagate the 
sample points along the discretized flow, and compute the expectation value at 
time t by averaging over ^E"* (a)-evaluations for the sample points. Our numerical 
experiments for Schrodinger equations in dimensions two and six, one with a 
torsional potential the other with a Henon-Heiles potential, indicate that both 
the initial sampling and the flow discretizations can acchieve sufficient accuracy 
to confirm the approximation's asymptotic convergence rate. 

1.1 Related research 

[LR10| has developed a discrete version of the Egorov theorem by propagating 
samples drawn from the initial Wigner function along the classical flow $*. 
The negativity of the Wigner function has been accounted for by stratified 
sampling as well as sampling from the uniform distribution on the effective 
support of the Wigner function. The observation, that the approximation order 
deteriorates from second to first order in e when merely replacing the Wigner 
by the Husimi function, has been formulated in |KLW091 Proposition 1] in the 
context of surface hopping algorithms. 

If the aim is not just the approximation of expectation values but of the 
wave function itself, then semiclassics have been successfully applied as well. 
The time-splitting algorithm of |FGL09] uses a Galerkin approximation with 
Hagedorn wave packets, which are products of polynomials with a Gaussian 
function of time- varying first and second moments. Related approximations are 
the Gaussian beam methods, whose convergence for the Schrodinger equation, 
the acoustic wave equation as well as strictly hyperbolic systems has been proven 
in }LRT10[ Theorem 1.1]. 

A complementary approximation ansatz combines Gaussians of fixed width 
with time-varying prefactors, whose evolution resembles the one of the width 
matrices of the Gaussian envelopes of both the Hagedorn wave packets and the 
Gaussian beams as well as the evolution of our correction term A* . The Herman- 
Kluk propagator SR09, Theorem 2], ^RlOi Theorem 1.2] as well as the frozen 
Gaussian approximation for strictly hyperbolic systems [LY12[ Theorem 4.1] 
have been systematically developed to fc-th order in e. 

1.2 Outline 

The following Section |2]begins with the comparison of Husimi and Wigner func- 
tions and provides asymptotic expansions in e for the transition from Weyl to 
Anti-Wick quantization and backwards in Lemma[l]and Lemma[2j respectively. 
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Section ^|3] develops and proves as Theorem [T] our first result, the second order 
Egorov approximation for Husimi functions. Our second result, Theorem [2] 
reformulates the correction term of Theorem [l] which has been derived as a 
time-integral along the <i>*-flow, with ordinary differential equations. Section ^ 
proposes a symmetric splitting method for the simultaneous discretization of 
A*, and r*. The final Section ^ presents our numerical results for two 
molecular Schrodinger equations in dimension two and six, respectively. 

2 Husimi functions and commutators 

Let ?/; : M'^ — > C be a square integrable function and W^{ip) : M^'' — > M, 



its Wigner function. The Wigner function is a continuous function on phase 
space. Its marginals are the position and momentum density of ?/;, respectively, 

W'mq,p)dp = / W{i^){q,p)dq = \T%ij){p)\\ 

where = (27re)~'^/^ J^^ e~*''^'/^'i/'('Z)'^'Z is the e-scaled Fourier trans- 

form. Moreover, associating with a Schwartz function a : M.'^'^ — )■ M the Weyl 
quantized operator 

(opW'=(a)V')(9) ^ (27^6)-" I I a{S±y,p)e^^^-yyp/^i^{y)dpdy, 



the Wigner function encodes corresponding expectation values via phase space 
integration, 

(^,opW<=(a)V)L2 = / a{z)W'mz)dz. 



However, despite its positive marginals, the Wigner function may attain negative 
values. This is easily seen for odd functions, which must satisfy >V'^(V')(0) = 



2.1 Husimi functions 

A proper smoothing of the Wigner function results in a nonnonegative phase 
space function, the so-called Husimi function: 

Definition 1. Let : M.'^ ^ C be square integrable and VV"^('0) its Wigner 
function. Let a : M^'' — > M &e a Schwartz function. We denote by 

Ge/2iz) = (^e)-'^e-l-l'/^ z € R'", 
the centered phase space Gaussian with covariance |Id. Then, 

are called the Husimi bmction ofip and the Anti-Wick operator of a, respectively. 
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By construction the Husimi function is a smooth function satisfying 



L2 



aiz)H'{^)iz) dz. 



(3) 



Especially, ~ \\'^\W2- One can prove, that H^itp) is the modu- 

lus squared of the FBI transform of ip and therefore nonnegative. Moreover, 
smoothing the Wigner function with Gaussians of smaller covariance does not 
yield positivity. In the following Proposition [T| we summarize these observa- 
tions combining arguments of the proofs of |Z08[ Proposition 2.20] and |B67[ 
Theorem 4.2]. 

Proposition 1. If <t > e, then yV^{ip) * > for all t/j G L'^{M.'^). Ifa<e, 
then there exists ijj G L^(IR'*) such that W^{ip) attains negative values. 

Proof. We write 

(W^(V')*G,/2) {q,p) 

= (7ra)-'^(27r£)-'^ 



r\y\'' 



dydx, 



smce 



We set a — {cr^/e^ — l) /(4tT). The change of variables x — = v, x + = w 
together with the relation 



and 
yields 

with 



\x~q\' + l\y\' = l\v^q\' + ^\w-q\' 



- + 4&|yP ^^\v~q\' + i,\w- qf + a\w - v\ 



{WW * G,/2) iq,p) = / f{v)fiw)e-''\^-''\"dvdw 



The Taylor series of a i-> e^°''""^ then gives 



li a — e, then a = and 



v^,■■■v,Jiv)e-''\'"^"dv 



f{v)dv 



> 0. 
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If CT > e, then G„i2 = * G(^„_^)/2 and 

If CT < e, then a < 0. Therefore, choosing '4>{v) = uie"'"'^, we obtain 

oo \k 



/c odd 



< 0. 



□ 



2.2 Weyl and Anti-Wick quantization 

The nonnegativity of the Husimi function is shared by the Anti-Wick quanti- 
zation, which satisfies a > op^^(a) > by identity (jsj). The foUowing 
Lemma expresses op'^^(a) as a Weyl quantized operator, whose symbol is ob- 
tained from a by adding higher order derivatives. It extends the first order 
approximation of |L101 Proposition 2.4.3] to higher orders. 

Lemma 1. Let a : M^'' — > M 6e a Schwartz function, n E N, and e > 0. There 
exists a Schwartz function rf^{a) : M.'^'^ — > M with sup^^Q ||op^°(r^(a))|| < oo 
such that 

where jk = j^r^fe)! f""^ k>0. 
Proof. We work with 

(a*G,/2)(^)-M-'^ 



a(C)e- 



-\C-AVe 



dc 



and Taylor expand a around the base point z, 

2n-l 

«(0 = E T-/'^{z){c-z,...x-z) 



k=0 



+ 



1 



(2n)! Jo 



(1 - (?)2"-ia(2")(z -I- 9{C - z)) dOiC - z, C - z). 



= :p(a(2"))(z,C) 

The terms with a-derivatives of odd degree do not contribute, since 
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for odd functions /. For the derivatives of even degree we obtain 
(tte)-'^ / a^^''\z){C - z, ...X - z)e-^'^''^'/' dC 

2d 

V "^^^ / y. ...y- e-\y\" dv 

2^ Ft.. L..^*i y*2.e ay 



-d k ^ d'^'^aiz) 



Since 



^1 5 ■ ■ •I'^k 



(tte)-'^ / C)(C - ^, C - z)e-l«"^l dC 



we obtain 



with 



n— 1 L. 



fe=0 



Moreover, 



||opW^«(a))|| < C sup 

iQi,i^j|<rd/2i+i 



for some constant C ~ C{d,n) > by the Calderon-Vaihancourt Theorem. □ 

Applying Lemma [l] repeatedly, we can also write op^°(a) as an Anti-Wick 
quantized operator by additively correcting a with its derivatives. 

Lemma 2. Let a : M^'* M. be a Schwartz function, n e N, and e > 0. There 
exists a Schwartz function Pn{a) : M?"^ — >■ M with sup^^Q ||op^''(p^(a))|| < oo 
such that 



opW<=(a) = opAW I „ _^ ^ e'^XkA'^a ] + £"opW<=(p^^(a)) 

with 



fc=i 



Xk 



l<m<A; 
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Proof. The first two appHcations of Lemma [T] yield 

op^%a) = op^'^ia) - opW° £"^7.1 A^^aj - e"opW^«(a)) 

(n—1 n—l — vi \ 
Vi — 1 V2 — I / 



vui = l 



The third apphcation of Lemma [T] results in 
op^-^Ca) 

V fe=i i;eN™.|ti|=fc y 

l<m<2 

E^' E nlr7..A'= 

/c=3 tieN3,|i,|=fc 

Cn-l 
k=i DeN'".|-u|=fe 

l<m<2 

By the (n — l)st application of Lemma [Tj we finally obtain 



/ 



op^°(a) = op 



where p^(a) 



AW 



E^' E (-irnr=i7..A^ 



T>2d 



fc=l i;eN™,|li|=fe 
l<m<fc 



is a Sehwartz function and 



+ £"opW°(p^,(a)), 



|op^'=(pf.(«))|| <C sup 

|a|,|/3|<rd/2l+l 



for some constant C = C{d,n) > 0. 



□ 



For further reference, we summarize the explicit form of the second order 
approximations of Lemma [l] and Lemma [2] 



op^W(a) ^ opW°(a+|Aa)+e2opW°(r|(a)), 
opW<=(a) = op^W(a- |Aa) +£2opWc(^e(^))_ 



(4) 
(5) 
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2.3 Commutators 

When approximating Husimi functions and Anti-Wick quantized expectation 
values for the solution of semiclassical Schrodinger equations, we will face com- 
mutators of Weyl and Anti-Wick quantized operators. For their asymptotic 
expansion we need additional notation: Let a,b E S{M.'^'^) be Schwartz func- 
tions. Then, the composition of their Weyl quantized operators is the Weyl 
quantized operator of the Moyal product ajjfo, 

opW°(a)opW<=(&) = opW<=(a^6), 

and the Moyal commutator ajJ6 — fefta has an e-expansion 

46 -b^a^Yl ^''^kia, b) + e"i?^(a, &), (6) 

A:=0 

where 

efe(a,6)(z) = ^ A{D)\a{z)b{z') - b{z)a{z'))\^,^^ (7) 
is obtained from the fc-fold application of 

A{D){a{z)b{z')) = i JVa(z) • V6(z'), J = 

and sup^>o l|oP^''('Rn('^' ^))ll < o^, see e.g. |Z121 Section 4.3]. We note, that 
Bfc(a, 6) = for fc even. For commutators of Weyl and Anti-Wick quantized 
operators we derive the following expansion. 

Lemma 3. Let a,b £ be Schwartz functions, n gN, and £ > 0. There 

exists a Schwartz function r^(a, b) : K^"* — >■ M with sup^^o W'^P^'^i'^'ni^^ ^))ll < 
such that 

n-l 

[opW-(a), opAW(6)] = eV^9k{a, b)) + £"opWe«(a, 6)) 

with 

Bk{a,b)=i -/iQmia,A^b), 

l+m=k+l 

where 0,„ and 7; have been defined in ^ and Lemma [1| respectively. 
Proof. By Lemma [T] 

J[opW<=(a),op^W(^)] 

n 

= JE^'t' [op^^(«),opW"(A'6)] +te" [opW°(a),opW<=(r^+^(6))] . 
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Using the Moyal expansion Q, we then obtain 
J[op^'=(a),op^W(6)] 

n n—l 

jEE^'^"7/op^'=(e„(a,A'6)) 

1=0 m=0 

/ n ^ 

n n—l 



1=0 m.=0 

FinaUy, we use Oo(a, &) = to rewrite 



n n — l 

;=0 m=0 k=0 l+m=k 

n-1 

= *E^'' E 7ze™(a,A'6). 

fc=0 i+m=fc+l 

□ 

Remark 1. We note that 9Q{a,b) — {a,b} and 9i{a,b) = ^{a,Ab}. Therefore, 
[opW<=(a),opAW(^)] ^ op'^^{{a,b+lAb})+eVirl{a,b)). 

3 Propagation 

In this section we apply the systematic transition between Weyl- and Anti-Wick 
calculus to prove a second order Egorov approximation for Anti-Wick quantized 
symbols. We can allow for unitary time evolutions e^'^*/"^ whose Hamiltonian 
H = op^°(/i) is obtained from a smooth symbol h : M^'* M of subquadratic 
growth. That is, for all I7I > 2 there exists > such that 

\\d-'h\\^ <C-,. 

In this situation, H is essentially self-adjoint in L^(]R'^) with core iS(R'^), 
^-iHt/e jg well-defined unitary operator on L^(K'^), and the Hamiltonian flow 

: M^*^ — >■ K^'^ associated with /i is a smooth mapping for all i G K. Schrodinger 
operators 

with smooth subquadratic potential V are one important example. 
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3.1 A first Egorov type approximation 

The proof of the Egorov type approximation is mainly buih on the following 
commutator estimate: 

Lemma 4. Let e > 0. Let b, c : M^'* — >■ M 6e a smooth function of subquadratic 
growth and a Schwartz function, respectively. Then, j [op^'^(6), op'^^(c)] is 
essentially self-adjoint in L^(M'^) with core 5(M'*), and there exists a Schwartz 
function c) : M^'' — > M with sup^^p lioP^'^(P2(^j c))|| < oo such that 

^[op^^(6),opAW(^)] 
e 

= op^"^ {{b- ^Ab,c} - ^tY{J D^bD^c)) + e^op'^%pl{b,c)). 
Proof. Since c is a Schwartz function, op^^(c) maps 

into and 

the commutator with op^°(6) is well-defined. By Lemma pi 

[op^^(6),op^W(c)] =opW'=({6,c+fAc})+eV<=(^^(6,c)) 
By Lemma [2] 

opW<=({6,c+fAc}) = op^'^{{b,c+lAc}^lA{b,c}) 

+ eV (p|({&,c+ |Ac}) - f;opAW(A{6, Ac}). 

By Lemma [l] 

op^W(A{6, Ac}) = opW°(A{6, Ac}) + eopW«(r^(A{6, Ac})), 
and therefore 



l [opW-(&), opAW(^)] ^ opAW(|^^ ^ _^ i^^i _ 1^1^^ ^ e'op'^^plib, c)) 

c) = r^(6, c) + p^b, c + f Ac}) - ^A{b, Ac} - ^rf (A{6, Ac}). 

a Schwartz function. We write A{6, c} = {Ab, c} + {b, Ac} + 2Y^^!Li{'^zkb, dz^c} 
and compute 

2d 2d 2d 



e 

with 



J2{dz,b,d,,c} = ^Ji?25(:,A:).Z?2c(:,fc)^^^2^(^^.)(^^2^)(.^^) 

fc=i 

= iv{JD%D^c), 



k=l fc=l fc=l 

-l2i, 7-i2 



where M(:, fc) and M{k, :) denote the fcth column and the fcth row of a matrix M, 
respectively. Therefore, 

{6,c+ |Ac}- |A{&,c} = {6- |A6,c}- f tr(Ji?26^2^)^ 
which concludes the proof. □ 
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Remark 2. If b is a polynomial of degree one, then the result of Lemma 
obviously simplifies to 

i[opW^(6),opAW(c)] =op^'^i{b,c}). 

If b is a polynomial of degree two, then the only simplification is due to {b — 
|A6, c} = {b,c}, since c) is not identical zero. 

Having prepared the necessary estimates, we state and prove our approxi- 
mation result for the unitary propagation of Anti-Wick operators. 

Theorem 1. Let h : M.'^'^ M. be a smooth function of subquadratic growth. 
Let a : W?'^ — > M &e a Schwartz function and t G M. There exists a constant 
C — C{a, h,t) > such that all e > 



with 



S*(a)=/" ti {J D^hD^ {a o <i>l)) o ^l-^ dT (8) 
Jo 

and $* : M^'* — > M^'' the Hamiltonian flow associated with h^ — h — |Aft,. 
Proof. We denote 

a(t) = ao$* - |S*(a) 

and observe that a{t) : M?'^ — > M is a Schwartz function. Since a(0) = a, we 
have 

^ [is (e''''^'0P^'^(°(^ - s))e-^"''') ds 
The second order commutator expansion of Lemma [4] yields 



[iI,op^W(a(t~.))] = {ov'^\h),op^'^{a{t~s))\ = 

op^w ^i^j^^ - s)} - I tT{JD^hD\{t - s))) + e^op^'^ip'^iKait - s))). 

Now, we compute 

dt{ao<i>l-^) = {h„ a o<S>l--} 

and 

dtE'-' (a) ^ tr {J D^hD^ao <Pl-^)) + {h„El'' (a)}. 

Therefore, 

9ta(t-s) = {h„a{t-s)}- |tr(jD2/iD2(ao$*-^)) 
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and 



^ [if,op^W(a(i - s))] - op^W(ata(i - s)) = 



The observation that 



ti{JD'^hD'^El-'{a)), pl{h,a{t- s)) eS{R^'^) 



concludes the proof. 



□ 



3.2 Corrections by ordineiry differential equations 

Our next task is to reformulate the time evolution of the correction term 



using ordinary differential equations which arc independent of the observable a. 

Lemma 5. Let a : M^'' — > M 6e a Schwartz function and h : M^'' — > M a smooth 
function of subquadratic growth. Let e > and : M^'' — )• R^"^ the Hamiltonian 
flow associated with kg = h — |A/i. Then, for all t e K, 



and r* : K^d ^ 

= / E {{JD'h)kidl,^^l)o^l--dr, 

where the Jacobian of is denoted as D^l = {V{^l)i, . . . ,^{^1)24)- 
Proof. We compute 





with A* : R^d ^ ^2dx2d ^ 




Then, 
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and 



_ »f 2d 



Adding the two terms gives the claimed identity. 



□ 



Having written the correction term S* (a) as the sum of a trace and an inner 
product, we can derive the desired a- independent ordinary differential equations 
for the time evolution of S* (a) . 

Theorem 2. Let a : M^'' R be a Schwartz function, h : M^'' — M a smooth 
function of subquadratic growth,, and t d R. Then, there exists a constant C = 
C{a, h,t) > such that for all e > 



with 



(a) = a o - I (tr(A* {D^a o $*)) + T* • (Va o $*)) 
where $* : M^'' -> M^'* is the Hamiltonian flow of h^ ^ h - |A/i, 

dfM ^JVheO^l, 

l^'^ solve 

dtAl = M,{t) + M,{t)Al+AlM,{tf, AO = 



and Ai 



^2d ^ |]j2dx2d pt 



p2d 



with 



dfTi 



M,{t) 



M,{t)Tl+tr{CdtyAl)t, 







(11) 

(12) 
(13) 



^2d 
T,2d 



^2dx2d 
j,2dx2d 



M,{t) = JD^ho^l, 
{C^{t)),u=dk{JD^h),jO<^\ 



Proof. We start with the reformulation of ^* (a) of Lemma [s] and observe that 
both relations ^ and (10) for A* and F*, respectively, involve time integrals of 
the form 



(/i(T)5/2(r))o$*--dr. 

"'0 

We compute 

dt f (/i(r)5/2(T)) o '^l-dT = /* dt ((/i(r)g/2(r)) o $^-) dr + /i(05/2(0 
Jo Jo 

and 

5t((/i(r)g/2(r))ocI>*-) 

= -9. ((/i(T)<?/2(r)) o $*--) + a, (/i(r).g/2(T)) o <i>*-^ 
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We obtain 



dt f\h{T)gh{T))o<i>l-dT 
Jo 



(/l(0)g/2(0))o$* + / dr{fl{T)gh{T))o<i>l-^dT. (14) 



This implies for 



Jo 



that 



rt 
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{{JD'^h, o <i>l){D%)'^ JD'^hD^l) o $*"^dT 

- I {{D^D^ JD^hD^liJD^h, o ^If) o $*"^dr, 
Jo 

= M,{t) + {JD^h, o $*)A* + K\{JD^h, o 

since the Hamihon equation (11) imphes drdz^^l — [j D^h^ o <I)^) dz^^l and 

= {JD^h,o^l){D'i>lf, 
drD^l = D<^>l{JD'^h,o^lf. 



Since J D^hf, o $* = + 0(e), the sohition A* of equation (12 1 satisfies 

A* = A* + 0(e), and we obtain 

a o $* - I S* (a) = a o - I (tr(A* [D^a o $*)) + P* • (Va o $*)) + 0{e^). 

For the time derivative of 

t 2d 
•^0 k,l=l 



equation (14) imphes 



We consider 9^^($pi = [(^J D'^h^ o $j) 9^^^$^) . and compute 



2d 



2d 

E iDK)in {DK)kj + {{JD^h, o • 
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Therefore, 



2d 



Jo 



j,k,l.n—l 



2d 



= {{JD^h,o<^l)f%+iT{C,{t)'^ll). 

Since J D^h^ o = M^{t) + 0(e), C,{t) = Ci{t) + 0(e), and A* = A* + 0(e), 
the solutfon r* of equation (13) satisfies T\—T\ + 0(e), and we conclude 



a o $* - I S* (a) = a o $* - I (tr(A* (i?^^ o $*)) + T* • (Va o $*)) + 0(e2). 

□ 

Corollary 1. Under the assumptions of Theorem^ there exists a constant 
C = C{a,h,t) > such that ipt = e~^^*'/^ipo satisfies 



Fl{a,){z)W{,l,o)(z) dz 



<Ce' 



(15) 



for all Tpa G L'^{M.'^) with H^oIIl^ — 1, where 

Fl{a) = ae o - I (tr(A* [D'^a o $*)) + T* • (Va o $*)) . 
Proof. We have 

(Vt,opW°(a)^t)i2 = (V,,op^^(a,)^t)L^ +0(e2) 
= (V'o,op^W(vI/*(ae))^o)L^+0(e2)= / vj,* (a,)(z) H^(V'o)(z) + 0(e2) 



and %{a^) = F*^ia) + 0{e^ 



□ 



4 Discretization 



For a numerical realization of the semiclassical approximation (151, we proceed 



in two steps. First, we discretize the phase space integral by equiweighted 
quadrature with nodes zi,...,Z]\[ £ R'^'^ distributed according to the Husimi 
function H^(iIjq). Then, we discretize the time evolution of $*, A*, and F* 
initialized with the sampling points zi , . . . , z^r . 
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4.1 Numerical quadrature 

To discretize the phase space integral 

we use, that the Husimi function of a normahzed wave function ipo is a proba- 
bihty density, and work either with Markov chain Monte Carlo quadrature or 
with Quasi-Monte Carlo quadrature. That is, we approximate 

1 ^ 

with quadrature nodes zi, . . . , zjv E M^"^ distributed according to ^{^{tpQ). 

4.1.1 Markov chain Monte Carlo quadrature 

We use Metropolis Monte Carlo and generate sample points zi , . . . , 2:jv , which 
form a Markov chain with stationary distribution W{tpQ), see e.g. |KLW09i 
§4.1]. If the chain is uniformly ergodic, then a central limit theorem holds, and 
there exists a constant 7 = 7(4'* (ce)) > such that for all c > 

lim F(\li^lia,))-Q^i^lia,))\ < ^) = ^ T e-^^'^dr. 
N^oo y vNJ v27r J_c 

4.1.2 Quasi-Monte Carlo quadrature 

For Gaussian wave packets ipo, the Husimi function 'H'^(^q) is a phase space 
Gaussian, and we view it as the density of a multivariate normal distribution. 
For notational simplicity, we focus on the special case of an isotropic Gaussian 
wave packet gz^ centered in zq — {qo,Po) € K^'', 

9z„iq) = (7r£)-'*/4exp (-ik - goP + ^Po ■ (q - go)) ■ (16) 

In this case, the Husimi function is the isotropic phase space Gaussian 

H%9z,){z) = (2^e)-'^ exp (-^|z - zoP) , 

that is, the density of a normal distribution with mean zq and covariance eld. 
Generating points of low star discrepancy with respect to the uniform distribu- 
tion on [0,1]^^^, that is, for example Sobol points, we use the cumulative dis- 
tributive function of the normal distribution and map them to points zi, . . . , zjy 
such that their star discrepancy with respect to the normal distribution 

2?*(zi, . . .,zn) 

= sup |^#{zj : Zj e (-oo,a),j = 1,...,A^} -'H''(-0o)((-oo,a))| 
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is optimal, that is, V* {zi, . . . , zn) = 0{{\og N)'^'^ /N). Then, the Koksma- 
Hlawka inequahty yields a constant 7 = j{'i>l{a^)) > such that 



\li^lia,))^QN{^liaM < 



(log NY" 

N '' 



where Cd > 2d dependends on the dimension of phase space, see |LR10[ 3.2]. 
4.2 A time splitting scheme 



For the simultaneous numerical solution of equations (111, (12), and (131, that 
is, for numerically computing $* , A* , and F* , we vectorize the Lyapunov matrix 
differential equation (12). For this purpose we use the Kronecker product 

aiiB ■ ■ ■ ainB 



of two matrices A e 



and B e 



a-nnB 



as well as the rowwise vectorization 



Vec(A) = (aii,ai2, . . . ,a„„) =Ae 



Lemma 6. The differential equations (11), (12), and (13) are equivalent to the 
coupled vector differential equation with Ad + Ad^ components 



dt I A* 




with 



^Id2d 
Ke{t) 
CS) MS) 



K,{t) = M,{t)®ld2d + lA2d®M,{t), 

■c,{tr- 



(17) 



Proof. First we show 



for A,B e M"^". 
Vec(A5)c„+i 



Vec(^B) ^{A® In) Yec{B) 
Indeed, for 1 < cn + i < 



acjbji — €S* In)(cn+i),{3n+i)Y'3c{B)jn+i 



^{A (g) In)(cn+^),k Vec(B)fe = {{A (g /„) Vec(B)) 



cn+2 
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Similarly one can prove that Yec{BA^) — {In <E) A)Vec{B). Hence, 
dtVec{Al) = M,{t)+Vec{M,{t)Al + AlM,{tf) 



Moreover, tr(Ci(t)A*) = C,{tYAl and 



□ 



For Schrodinger Hamiltonians h{q,p) = + V{q) with smooth sub- 



quadratic potential V -.R —, 

MS) = JD^ho^l 



we observe that 




Id. 



while Ci{t)jk — dk{JD^h)ij o satisfies Ci{t) = for i = 1, . . . , d and depends 
on third derivatives of V and {^l)q ioi i = d+1, . . . , 2d. This motivates to split 





1 ' \ 









A* 


- A^ ' 





\n J 


V r* J 


V / 



and to rewrite the differential equation (17) as 



5tT* = A{Tl)T\ + iV(T^) + BT{ 



(18) 



with 











' ) 

























, B = 




















Keit) 






















Ceit) 


Me{t)J 












0/ 



and 



7V(T*) 



v 





-VF(($*),) 
M,{t) 




Let 0^ and (j>l be the flows of the systems 



^tT* = A(T*)T* +Ar(T*), dtT^^BTl, 



19 



respectively, and 

the Strang splitting for /i > 0, which provides a second order time discretization 



of equation 18 see e.g. |HWL06[ Section III. 3.4]. This splitting scheme produces 



an approximate solution via 

= T° + |(yl(TO)T; + iV(T°)), 

= Ti^V |(A(T^)Ti/^ + A^(T^)). 

and discretizes $* by the Stormer-Verlet scheme, while A* and F* are discretized 
by the midpoint rule. 

Remark 3. The flows and <j)\ are numerically easy to evaluate, since they 
are determined by differential equations with constant coefficient matrices and 
a constant inhomogeneity. 



5 Numerical experiments 

Let : M'' ^ C be the solution of the Schrodinger equation ^ and a : M^'^ R 
a Schwartz function. The preceeding algorithmical considerations suggest 

1 ^ ~ 

(V'*,opW''(a)V't)^. « -5]i^*(a)(z,), (19) 

where Zi, . . . ,Zn S M^'^ are distributed according to 'H'^{tpo)- F*{a) is obtained 
by lt/h\ iterates of the Strang splitting scheme F'^ and the application of the 
resulting vector to a and its derivatives according to Corollary [l] Our numerical 
experiment^validate this approach for different observables a: the position and 
momentum operators defined by 

as well as the potential, kinetic, and total energy operators defined by 

{q,p)^V{q), iq,p)^^\p\^, {q,p) ^ h{q,p), 

respectively. The main focus of the experiments is on the positive answer of the 
following two questions: Is the algorithm feasible in a moderate high dimensional 
setting? Can the asymptotic O(e^) accuracy be acchieved efficiently, that is, 
with a reasonable number of sampling points N and time stepping h? 

^AU experiments have been performed with MATLAB 7.14 on a 3.33 GHz Intel Xeon 
X5680 processor. 
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5.1 Six dimensions 



Our first experiment is concerned with the Henon-Heiles potential in six dimen- 
sions, 

6 5 

with cr* l/VSO. As in |FGL09[ §5.4] and |LR10[ §6], we choose the semiclassi- 
cal parameter e — 0.01 and a Gaussian initial state ( [l6| centered in zq — (qo, 0) 
with qo = (2,2,2,2,2,2)""". Since a grid based reference solution of the six- 
dimensional Schrodinger equation is not feasible, we compare the following three 
asymptotic particle methods: 



The second order Husimi based method defined in ( 19 1 using N transformed 
Sobol points for the numerical quadrature. 



B. The naive Husimi method 

1 ^ 

(Vt,opW°(a)^,)^,«-^(aoci>*)(z,-), 



where <I>* is the Stormer-Verlet discretization of q = p, p = —W{q). The 
quadrature nodes zi , . . . , z^v are obtained by transforming Sobol points, 
such that they are distributed according to 7i^{ipo). This approach is first 
order accurate with respect to e. 

C. The second order Wigner based method 

1 ^ 

(^„opW^(a)V;,)^.«-5](ao$*)(^^.) 

of |LR101 §6], where the quadrature nodes wi, . . . ,wn are obtained by 
transforming Sobol points, such that they are distributed according to the 
initial Wigner function yV"^(?/'o)- 

For the six-dimensional experiments in this paragraph we use N = 2 ■ 10** 
sampling points and the time stepping h = 5- 10~^ to reach an accuracy confirm- 
ing our theoretical expectations. Figure [T] illustrates that algorithms A and C 
are second order accurate with respect to e, while the naive Husimi approach 
of algorithm B is only first order with respect to e. The computing time for 
the algorithms B and C, which use uncorrected classical transport, is 75 sec- 
onds, while the corrected Husimi algorithm A requires 140 minutes due to the 
evaluations of the second and third derivatives of the potential V as well as the 
handling of the vector ($*, A*,r*), which has length Ad + AcP = 168. 
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Figure 1: The difference of the expectation values of the kinetic, potential 
and total energy for the six-dimensional Henon-Heiles system computed by al- 
gorithms A and C (left) as well as B and C (right) as functions of time. The 
plots confirm that the corrected Husimi method A is second order accurate with 
respect to £, while the naive Husimi approach B is only first order accurate. 



5.2 Different semiclassical parameters 

The next experiments are performed for a two-dimensional torsional potential, 

F(9) = 2-cos(gi)-cos(g2), (20) 

and Gaussian initial data g^g centered in zq = (qoiOiO)"^, while varying the 
semiclassical parameter e. 



e = 10 ^. The position center is — (0, 1). The discretization (191 uses N 
10"^ Sobol points and the time stepping h = lO^^. 



e = 10 2. The position center is qo — (0.5,0.8). The discretization (19) uses 
= 3 • 10"* Sobol points and the time stepping h = 10^^. 



— in-3 



The position center is qo ~ (1,0). The discretization (19 1 uses A^ 



10- 

8 • 10^ Sobol points and the time stepping h = 5 ■ 10" 



e 


#tiniesteps 


domain 




space grid 


10-' 


5 • 103 


[-2,2] X [-2, 


2] 


2048 X 2048 


10-2 


5 • 103 


[-2,2] X [-2, 


2] 


2048 X 2048 


10-3 


104 


[-2,2] X [-1, 


1] 


4096 X 2048 



Table 1: The discretization parameters for the grid-based reference solutions 
computed by a Strang splitting scheme with Fourier collocation. 



We generate numerically converged reference solutions for the Schrodinger 
equation ([T]) using a Strang splitting scheme with Fourier collocation for the 
discretization of the Laplacian. The discretization parameters for the reference 
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solutions, that is, the number of time steps, the computational domain, and the 
size of the space grid, are summarized in Tablejl] The comparison of the Husimi 
algorithm [19] with the expectation values inferred from the reference solutions 
is presented in Figure [2j It confirms our expectations and shows second order 
accuracy with respect to e. 




<i ' ' ' ' 1 

4 8 12 16 20 



position 

_ momentum 

kinetic 
potentiai 
total 











4 8 12 16 20 

time 



Figure 2: The errors of the expectation values of position, momentum, po- 
tential, kinetic, and total energy for the two-dimensional torsional potential 
computed with the corrected Husimi algorithm ( [l9| as functions of time. The 
semiclassical parameter is chosen as e = 10~^ (top), e — 10~^ (middle), and 
e — 10^^ (bottom). In all three cases, the errors are of the order 
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5.3 Superposition of Gaussian wave packets 

The last set of experiments studies the dynamics of superpositions of Gaus- 
sian wave packets for the Schrodinger equation ([T]) with the torsional potential 



defined in (20). The initial data ipo are superpositions + with phase 
space centers zi,Z2 G M^**, d — 2, normalized such that HV'oIIl^ — 1- Since 
the torsional potential is symmetric around the origin, we consider a symmetric 
placement of the centers zi and Z2 as well as an asymmetric one. For Gaussian 
superpositions the Husimi function can easily be calculated as 

H'i9z,+gz,)=H%gz,) + H%gz,) + 2Cz„z,, 
where the cross term 

Cz,,z,{z) = (27r£)-'*exp {-j^\z_\^) exp {-^\z - z+p) cos (^(ci,2 - Jz ■ z_)) 

has a Gaussian envelope centered around the mean — \{zi+Z2) and oscillates 
with a frequency proportional to the difference z_ — zi — Z2. The shift is 
Ci.2 — <l{zi)p{zi) — q{z2)p{z2)- The cross term contains a constant damping 
factor, which is exponentially small in |z_p. This allows the tails of 'H^(gzi) 
and H'^ (g^^ ) to absorb the oscillations of the cross term, which turns the Husimi 
function positive. Our experiments are the following: 

£ — 10^^. The first center is zi = (1,0,0,0)"^, while the second center is Z2 = 
— zi or Z2 = (—0.4, —0.2, 0, 0)'^. In both cases, we generate Markov chains 
by an Metropolis algorithm with jumps between sampling regions centered 
around Zi, Z2, and z+, see |KLW091 §4.1]. The symmetric superposition 
is sampled with = 3 • 10^ points, while the asymmetric case requires 
N = 5 ■ 10^ sampling points. The final results are arithmetic means over 
ten independent runs of this set up, which provides unbiased estimates for 
the phase space integrals. 

e = 10^^. The first center is zi = (1, 0.5, 0, 0)"^, while the second center is Z2 = 
— Zi or Z2 = (-0.4,-0.45,0,0)"'". In both cases, the cross term of the 
Husimi function is smaller than 10^^^ due to the exponentially damping 
factor exp(— g^|z_p). Therefore, the cross term is neglected, and the two 
Gaussians H'^lgzi) and H'^igzi) si's each sampled by = 5 • 10^ Sobol 
points. 

All experiments use the time stepping h — 10^^ for the corrected flow F'\ 
The comparison with the grid based reference solutions obtained with the dis- 
cretization parameters of Table [T] is presented in Figure |3] Again, it confirms 
second order accuracy with respect to the semiclassical parameter. We mention, 
that for the symmetrically placed phase space centers zi and Z2, the position 
and momentum expectations are more accurately approximated than the energy 
expectations. 
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Figure 3: The errors of expectation values for the two-dimensional torsional 
potential with e = 10~^ (top row) and e = 10~^ (bottom row) computed with 
the corrected Husimi algorithm ( |l9| ). The initial data are superpositions of 
Gaussians wave packets placed symmetrically (left column) or asymmetrically 
(right column). In all four cases, the errors are of the order e^. 



A Appendix: Composition formula 

The calculus developed in § [2] implies a composition formula for Anti-Wick 
quantized operators. 

Lemma 7. Let a, b : M^^ R be Schwartz functions, n E N, and e > 0. Then, 
opAW(a)opAW(6) = op^w (ab+ e™c„J + 0{e"+') 
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with 



\i=l / k+l+j=s 



\v\+s—m,vGN^ ,\v\<n — s 
or v—0 and r— 



/or m > 1, where Qj{a, b) is defined 
Proof. We use Lemma [T] and write 



Vfe=0 



Now we apply Q and obtain 
op^W(a)opAW(^) 



- op^M E e'+'(7fcA'=a)tt(7iA'6) +0(e"+i) 



op 



[ E e'+'+-'ej(7;cAV7/A'6) I 

\k+l+j<n J 



= Ee^op"^' E e,(7.AV7/A'fe) 

m=o yfc+i+j=m y 

We set X]fc+;+j=m 0j(7feA'''a, 7iA'5) =: dm and observe that by Lemma[2] 
opW^ (dm) 



= op 



,AW 



V 



vGN'^ ,\v\<n—m 
l<r<n—m 



Collecting all the terms of the same order in e yields 
with 



m— 1 



E (-ir n^- ^""'^^ 



or v—0 and r— 



Remark, for n — 1, we obtain 

op^W(a)op^W(^) ^ opAW + i| JVa • V6 - fVa • V6) + ©(e^) 
whose real part coincides with ILIO, Lemma 2.4-6]. 
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